با استفاده از داده های زلزله ها در ایران و جهان به سوالات زیر پاسخ دهید.
library(readr)
library(dplyr)##
## Attaching package: 'dplyr'## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionlibrary(highcharter)## Highcharts (www.highcharts.com) is a Highsoft software product which is## not free for commercial and Governmental uselibrary(ggplot2)## Warning: package 'ggplot2' was built under R version 3.4.4library(plotly)## Warning: package 'plotly' was built under R version 3.4.4##
## Attaching package: 'plotly'## The following object is masked from 'package:ggplot2':
##
## last_plot## The following object is masked from 'package:stats':
##
## filter## The following object is masked from 'package:graphics':
##
## layoutlibrary(gganimate)
library(ggmap)## Google Maps API Terms of Service: http://developers.google.com/maps/terms.## Please cite ggmap if you use it: see citation("ggmap") for details.##
## Attaching package: 'ggmap'## The following object is masked from 'package:plotly':
##
## windlibrary(h2o)## Warning: package 'h2o' was built under R version 3.4.4##
## ----------------------------------------------------------------------
##
## Your next step is to start H2O:
## > h2o.init()
##
## For H2O package documentation, ask for help:
## > ??h2o
##
## After starting H2O, you can use the Web UI at http://localhost:54321
## For more information visit http://docs.h2o.ai
##
## ----------------------------------------------------------------------##
## Attaching package: 'h2o'## The following objects are masked from 'package:stats':
##
## cor, sd, var## The following objects are masked from 'package:base':
##
## %*%, %in%, &&, ||, apply, as.factor, as.numeric, colnames,
## colnames<-, ifelse, is.character, is.factor, is.numeric, log,
## log10, log1p, log2, round, signif, trunclibrary(rworldmap)## Warning: package 'rworldmap' was built under R version 3.4.4## Loading required package: sp## Warning: package 'sp' was built under R version 3.4.4## ### Welcome to rworldmap ##### For a short introduction type : vignette('rworldmap')library(stringr)## Warning: package 'stringr' was built under R version 3.4.4equake = read_csv("E:\\edu\\TahlilDade\\data\\worldwide.csv")## Parsed with column specification:
## cols(
## .default = col_double(),
## time = col_datetime(format = ""),
## magType = col_character(),
## nst = col_integer(),
## net = col_character(),
## id = col_character(),
## updated = col_datetime(format = ""),
## place = col_character(),
## type = col_character(),
## magNst = col_integer(),
## status = col_character(),
## locationSource = col_character(),
## magSource = col_character()
## )## See spec(...) for full column specifications.iequake = read_rds("E:\\edu\\TahlilDade\\data\\iran_earthquake.rds")
disaster = read_delim("E:\\edu\\TahlilDade\\data\\disaster.txt", "\t", escape_double = FALSE, trim_ws = TRUE)## Parsed with column specification:
## cols(
## .default = col_integer(),
## FLAG_TSUNAMI = col_character(),
## SECOND = col_character(),
## EQ_PRIMARY = col_double(),
## EQ_MAG_MW = col_double(),
## EQ_MAG_MS = col_double(),
## EQ_MAG_MB = col_character(),
## EQ_MAG_ML = col_double(),
## EQ_MAG_MFA = col_character(),
## EQ_MAG_UNK = col_double(),
## COUNTRY = col_character(),
## STATE = col_character(),
## LOCATION_NAME = col_character(),
## LATITUDE = col_double(),
## LONGITUDE = col_double(),
## MISSING = col_character(),
## DAMAGE_MILLIONS_DOLLARS = col_character(),
## TOTAL_MISSING = col_character(),
## TOTAL_MISSING_DESCRIPTION = col_character(),
## TOTAL_DAMAGE_MILLIONS_DOLLARS = col_character()
## )
## See spec(...) for full column specifications.sequake = read_rds("E:\\edu\\TahlilDade\\data\\historical_web_data_26112015.rds")۱. با استفاده از داده های historical_web_data_26112015.rds و استفاده از نمودار پراکنش سه بعدی بسته plotly نمودار طول، عرض و عمق زلزله ها را رسم نمایید. علاوه بر آن بزرگی هر نقطه را برابر بزرگی زمین لرزه قرار دهید.
plot_ly(sequake, x=~Longitude, y=~Latitude, size=~Magnitude, z=~Depth, sizes=c(1,1000))## No trace type specified:
## Based on info supplied, a 'scatter3d' trace seems appropriate.
## Read more about this trace type -> https://plot.ly/r/reference/#scatter3d## No scatter3d mode specifed:
## Setting the mode to markers
## Read more about this attribute -> https://plot.ly/r/reference/#scatter-mode## Warning: package 'bindrcpp' was built under R version 3.4.4۲. پویانمایی سونامی های تاریخی را بر حسب شدت بر روی نقشه زمین رسم نمایید.(از داده زلزله های بزرگ استفاده نمایید.)
bbox <- c(left = -170, bottom = -60, right = 170, top = 80)
World_map <- get_stamenmap(bbox, zoom = 3, maptype="terrain")## 48 tiles needed, this may take a while (try a smaller zoom).## Source : http://tile.stamen.com/terrain/3/0/0.png## Source : http://tile.stamen.com/terrain/3/1/0.png## Source : http://tile.stamen.com/terrain/3/2/0.png## Source : http://tile.stamen.com/terrain/3/3/0.png## Source : http://tile.stamen.com/terrain/3/4/0.png## Source : http://tile.stamen.com/terrain/3/5/0.png## Source : http://tile.stamen.com/terrain/3/6/0.png## Source : http://tile.stamen.com/terrain/3/7/0.png## Source : http://tile.stamen.com/terrain/3/0/1.png## Source : http://tile.stamen.com/terrain/3/1/1.png## Source : http://tile.stamen.com/terrain/3/2/1.png## Source : http://tile.stamen.com/terrain/3/3/1.png## Source : http://tile.stamen.com/terrain/3/4/1.png## Source : http://tile.stamen.com/terrain/3/5/1.png## Source : http://tile.stamen.com/terrain/3/6/1.png## Source : http://tile.stamen.com/terrain/3/7/1.png## Source : http://tile.stamen.com/terrain/3/0/2.png## Source : http://tile.stamen.com/terrain/3/1/2.png## Source : http://tile.stamen.com/terrain/3/2/2.png## Source : http://tile.stamen.com/terrain/3/3/2.png## Source : http://tile.stamen.com/terrain/3/4/2.png## Source : http://tile.stamen.com/terrain/3/5/2.png## Source : http://tile.stamen.com/terrain/3/6/2.png## Source : http://tile.stamen.com/terrain/3/7/2.png## Source : http://tile.stamen.com/terrain/3/0/3.png## Source : http://tile.stamen.com/terrain/3/1/3.png## Source : http://tile.stamen.com/terrain/3/2/3.png## Source : http://tile.stamen.com/terrain/3/3/3.png## Source : http://tile.stamen.com/terrain/3/4/3.png## Source : http://tile.stamen.com/terrain/3/5/3.png## Source : http://tile.stamen.com/terrain/3/6/3.png## Source : http://tile.stamen.com/terrain/3/7/3.png## Source : http://tile.stamen.com/terrain/3/0/4.png## Source : http://tile.stamen.com/terrain/3/1/4.png## Source : http://tile.stamen.com/terrain/3/2/4.png## Source : http://tile.stamen.com/terrain/3/3/4.png## Source : http://tile.stamen.com/terrain/3/4/4.png## Source : http://tile.stamen.com/terrain/3/5/4.png## Source : http://tile.stamen.com/terrain/3/6/4.png## Source : http://tile.stamen.com/terrain/3/7/4.png## Source : http://tile.stamen.com/terrain/3/0/5.png## Source : http://tile.stamen.com/terrain/3/1/5.png## Source : http://tile.stamen.com/terrain/3/2/5.png## Source : http://tile.stamen.com/terrain/3/3/5.png## Source : http://tile.stamen.com/terrain/3/4/5.png## Source : http://tile.stamen.com/terrain/3/5/5.png## Source : http://tile.stamen.com/terrain/3/6/5.png## Source : http://tile.stamen.com/terrain/3/7/5.pngdisaster %>% filter(FLAG_TSUNAMI=='Tsu') -> Tsu_disaster
ggmap(World_map)+
geom_point(aes(x=LONGITUDE, y=LATITUDE), data = Tsu_disaster)+
transition_states(
EQ_PRIMARY,
transition_length = 2,
state_length = 1
) +
enter_fade() +
exit_fade()## Warning in (function (d, t, en, ex, es) : NAs introduced by coercion## Warning: Removed 79 rows containing missing values (geom_point).
## Warning: Removed 79 rows containing missing values (geom_point).
۳. نمودار چگالی دو بعدی زلزله های تاریخی ایران را رسم کنید.( از داده iran_earthquake.rds و لایه stat_density_2d استفاده نمایید).
iran = read_rds("E:\\edu\\TahlilDade\\data\\iran_earthquake.rds")
iran_map = read_rds("E:\\edu\\TahlilDade\\data\\Tehrn_map_6.rds")
ggmap(iran_map) + stat_density_2d(data=iran, aes(x=Long, y=Lat))## Warning: Removed 243 rows containing non-finite values (stat_density2d).
۴. احتمال اینکه در ایران در پنج سال آینده زلزله به بزرگی هفت ریشتر رخ دهد را محاسبه کنید. (از احتمال شرطی استفاده کنید.)
iequake$year = iequake$OriginTime %>% str_extract("(\\d)*")
iequake %>%
filter(Mag > 7) %>%
.$year %>%
unique() %>%
length() -> year_with_quake
p = 1 - year_with_quake/5
print(1 - p^5)## [1] 0.98976۵. بر اساس داده های زلزله های بزرگ ابتدا تعداد و متوسط کشته زلزله ها را بر حسب کشور استخراج نمایید. سپس نمودار گرمایی تعداد کشته ها را بر روی کره زمین رسم نمایید.(مانند مثال زیر!)
disaster %>%
group_by(COUNTRY) %>%
summarise(sum_death = sum(TOTAL_DEATHS, na.rm = TRUE),
mean_death = mean(TOTAL_DEATHS, na.rm = TRUE)) -> death_quake
matched <- joinCountryData2Map(death_quake, joinCode="NAME", nameJoinColumn="COUNTRY")## 138 codes from your data successfully matched countries in the map
## 16 codes from your data failed to match with a country code in the map
## 105 codes from the map weren't represented in your datamapCountryData(matched, nameColumnToPlot="sum_death")## You asked for 7 quantiles, only 4 could be created in quantiles classification
۶. با استفاده از داده لرزه های بزرگ و به وسیله طول، عرض، شدت، عمق مدلی برای پیش بینی تعداد کشته های زلزله بیابید.
h2o.init()## Connection successful!
##
## R is connected to the H2O cluster:
## H2O cluster uptime: 1 days 4 hours
## H2O cluster version: 3.16.0.2
## H2O cluster version age: 7 months and 20 days !!!
## H2O cluster name: H2O_started_from_R_ASUS_jrj465
## H2O cluster total nodes: 1
## H2O cluster total memory: 3.34 GB
## H2O cluster total cores: 8
## H2O cluster allowed cores: 8
## H2O cluster healthy: TRUE
## H2O Connection ip: localhost
## H2O Connection port: 54321
## H2O Connection proxy: NA
## H2O Internal Security: FALSE
## H2O API Extensions: Algos, AutoML, Core V3, Core V4
## R Version: R version 3.4.3 (2017-11-30)## Warning in h2o.clusterInfo():
## Your H2O cluster version is too old (7 months and 20 days)!
## Please download and install the latest version from http://h2o.ai/download/disaster %>%
select(LATITUDE, LONGITUDE, FOCAL_DEPTH, EQ_PRIMARY, DEATHS) %>%
mutate(DEATHS = ifelse(!is.na(DEATHS),DEATHS,0)) %>%
.[complete.cases(.),] %>%
as.h2o()-> train##
|
| | 0%
|
|=================================================================| 100%h2o.glm(x=colnames(train), y=c('DEATHS'), train)## Warning in .verify_dataxy(training_frame, x, y): removing response variable
## from the explanatory variables##
|
| | 0%
|
|=================================================================| 100%## Model Details:
## ==============
##
## H2ORegressionModel: glm
## Model ID: GLM_model_R_1532021741235_2
## GLM Model: summary
## family link regularization
## 1 gaussian identity Elastic Net (alpha = 0.5, lambda = 1.5598 )
## number_of_predictors_total number_of_active_predictors
## 1 4 4
## number_of_iterations training_frame
## 1 1 file6053bb587_sid_966d_2
##
## Coefficients: glm coefficients
## names coefficients standardized_coefficients
## 1 Intercept -2128.626172 872.312939
## 2 LATITUDE 9.903106 248.426544
## 3 LONGITUDE -0.395260 -36.942863
## 4 FOCAL_DEPTH -2.618074 -183.862533
## 5 EQ_PRIMARY 457.651819 501.255831
##
## H2ORegressionMetrics: glm
## ** Reported on training data. **
##
## MSE: 97262470
## RMSE: 9862.174
## MAE: 1612.495
## RMSLE: NaN
## Mean Residual Deviance : 97262470
## R^2 : 0.008045576
## Null Deviance :293271587761
## Null D.o.F. :2990
## Residual Deviance :2.90912e+11
## Residual D.o.F. :2986
## AIC :63513.33۷. با استفاده از داده worldwide.csv به چند سوال زیر پاسخ دهید. تحقیق کنید آیا می توان از پیش لرزه، زلزله اصلی را پیش بینی کرد؟
equake$place %>%
str_locate('of') %>%
.[,2] %>%
ifelse(is.na(.), -1, .) -> indexes
indexes <- indexes + 2
equake$place1 = indexes %>% str_sub(equake$place, start = .)
equake %>%
group_by(place1) %>%
arrange(time) %>%
mutate(nex = lead(mag)) %>%
filter(nex > mag) %>%
lm(nex~mag, .) %>%
summary()##
## Call:
## lm(formula = nex ~ mag, data = .)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.4878 -0.2907 -0.1150 0.1265 4.9022
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.619510 0.011730 52.81 <2e-16 ***
## mag 0.951314 0.003178 299.31 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.408 on 25182 degrees of freedom
## Multiple R-squared: 0.7806, Adjusted R-squared: 0.7806
## F-statistic: 8.958e+04 on 1 and 25182 DF, p-value: < 2.2e-16۸. گزاره " آیا شدت زلزله به عمق آن بستگی دارد" را تحقیق کنید؟ (طبیعتا از آزمون فرض باید استفاده کنید.)
با استفاده از ANOVA بله
equake %>% mutate(new_mag = as.integer(mag)) %>%
aov(depth~new_mag, .) -> a
summary.aov(a)## Df Sum Sq Mean Sq F value Pr(>F)
## new_mag 1 37581665 37581665 3136 <2e-16 ***
## Residuals 60629 726613331 11985
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1۹. میانگین سالانه زلزله ها را بر حسب کشور به دست آورید. آیا میتوان دلیلی در تایید یا رد تئوری هارپ ارائه کرد.
متاسفانه این تیوری را نمی توان رد یا اثبات کرد روال کار این بود که سعی کنیم ببینیم آیا میانگین تعداد زلزله های کشور ها نسبت به سال فرقی می کند یا نه
equake$country = equake$place %>% str_extract(", (\\w)*") %>% str_sub(start = 3L)
equake$year = equake$time %>% str_extract("(\\d)*")
equake %>%
group_by(country, year) %>%
summarise(meanmag = mean(mag, na.rm = TRUE)) -> mean_quake
aov(year~meanmag,mean_quake) -> a
summary.aov(a)## Df Sum Sq Mean Sq F value Pr(>F)
## meanmag 1 0.5 0.4535 0.431 0.512
## Residuals 574 604.5 1.0532۱۰. سه حقیقت جالب در مورد زلزله بیابید.
۱- در ایران از ژاپن بیشتر زلزله می آید!
disaster %>%
filter(YEAR > 2000) %>%
group_by(COUNTRY) %>%
summarise(num = n()) %>%
filter(COUNTRY %in% c('IRAN', 'JAPAN')) %>% View۲- شدت زلزله به سال آن بستگی دارد! در واقع این احتمالا ناشی از تغییر توانایی بشر در ثبت زلزله هاست
disaster %>%
mutate(DES_YEAR = as.factor(as.integer(YEAR/100))) %>%
aov(EQ_PRIMARY~DES_YEAR, .) %>% summary.aov()## Df Sum Sq Mean Sq F value Pr(>F)
## DES_YEAR 28 428 15.3 15.3 <2e-16 ***
## Residuals 4208 4208 1.0
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 1789 observations deleted due to missingness۳- باوجود اینکه تعداد کشته ها باید به ساعت زلزله بستگی شدیدی داشته باشد ولی این موضوع را نمی توان ثابت کرد.
disaster %>%
aov(TOTAL_DEATHS~HOUR, .) %>% summary.aov()## Df Sum Sq Mean Sq F value Pr(>F)
## HOUR 1 3.686e+07 36861544 0.135 0.714
## Residuals 1263 3.453e+11 273422883
## 4761 observations deleted due to missingness